% Consider any of the estimated models and run this script as a supplement

% Re-estimating P dynamics
[theta2Step3,AVarTheta2Step3] =step2SR_threshold_macro(AVarTheta11Step3,outputStep3,bandwidthTstep2,setupStep3.econAct,h0RegimeOn,hxRegimeOn,numBootStep2);

% Enforce consistency of sigma between Q and P
posTheta12Step2 = zeros(1,length(namesTheta12));
for i=1:length(namesTheta12)
    theta2Step3.(namesTheta12{i}) = theta2Step2.(namesTheta12{i});
    posTheta12Step2(1,i) = find(strcmp(AVarTheta2Step3.names,namesTheta12(i)));
end
AVarTheta2Step3.VarRobust(posTheta12Step2,posTheta12Step2) = AVarTheta2Step2.VarRobust(posTheta12Step2,posTheta12Step2); 

% Mimicking a two-state Markov-model
% Expansions
theta2Step3.h03_1  = mean(outputStep3.factors(1,setupStep3.econAct>=0),2);
theta2Step3.h04_1  = mean(outputStep3.factors(2,setupStep3.econAct>=0),2);
theta2Step3.hx31_1 = 0;
theta2Step3.hx32_1 = 0;
theta2Step3.hx33_1 = 0;
theta2Step3.hx34_1 = 0;
theta2Step3.hx41_1 = 0;
theta2Step3.hx42_1 = 0;
theta2Step3.hx43_1 = 0;
theta2Step3.hx44_1 = 0;

% Recessions
theta2Step3.h03_2  = mean(outputStep3.factors(1,setupStep3.econAct<0),2);
theta2Step3.h04_2  = mean(outputStep3.factors(2,setupStep3.econAct<0),2);
theta2Step3.hx31_2 = 0;
theta2Step3.hx32_2 = 0;
theta2Step3.hx33_2 = 0;
theta2Step3.hx34_2 = 0;
theta2Step3.hx41_2 = 0;
theta2Step3.hx42_2 = 0;
theta2Step3.hx43_2 = 0;
theta2Step3.hx44_2 = 0;

theta2Step3.sigma31 = 0;
theta2Step3.sigma32 = 0;
theta2Step3.sigma33 = 0;

theta2Step3.sigma41 = 0;
theta2Step3.sigma42 = 0;
theta2Step3.sigma43 = 0;
theta2Step3.sigma44 = 0;


% Turning off feedback from yield curve to PMI
%theta2Step3.gamax1 = 0;
%theta2Step3.gamax2 = 0;
%theta2Step3.gamax3 = 0;
%theta2Step3.gamax4 = 0;

% Simulating moments
momentsData             = getMomentsData(loadData,outputStep3,thresholdLevel);
momentsModel            = simulator_threshold(theta2Step3,outputStep3,thresholdLevel,numSim);

% Plotting
plotUnconditionalMoments(momentsModel,momentsData)
plotForecastRegressions(momentsModel,momentsData)

% Fit of bond yields in the data
[T,ny] = size(setupStep3.yields);
yieldsModel = zeros(T,ny);
model = outputStep3.model;
model.factorSignResOn = 0;
for t=1:T
    if setupStep3.econAct(t,1) >= 0
        x2 = [theta2Step3.h03_1;theta2Step3.h04_1];
    else
        x2 = [theta2Step3.h03_2;theta2Step3.h04_2];
    end
    fit = modelFit(model,x2,t);
    yieldsModel(t,:) = fit';
end
residuals = setupStep3.yields-yieldsModel;
RMSEresiduals = sqrt(mean(residuals.^2,1));
plot([yieldsModel(:,39),setupStep3.yields(:,39)])

filename = 'MarkovSwitchingProxy';
 save([pwd,'\Results\',filename,'.mat'],...
     'theta11Step3','outputStep3','setupStep3','AVarTheta11Step3','theta2Step3',...
     'momentsModel','momentsData')
